Kinetic roughening with anisotropic growth rules 
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Inspired by the chemical etching processes, where experiments show that growth rates depending 
on the local environment might play a fundamental role in determining the properties of the etched 
surfaces, we study here a model for kinetic roughening which includes explicitly an anisotropic effect 
in the growth rules. Our model introduces a dependence of the growth rules on the local environment 
conditions, i.e. on the local curvature of the surface. Variables with different local curvatures of 
the surface, in fact, present different quenched disorder and a parameter p (which could represent 
different experimental conditions) is introduced to account for different time scales for the different 
classes of variables. We show that the introduction of this time scale separation in the model leads 
to a cross-over effect on the roughness properties. This effect could explain the scattering in the 
experimental measurements available in the literature. The interplay between anisotropy and the 
cross-over effect and the dependence of critical properties on parameter p is investigated as well as 
the relationship with the known universality classes. 
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I. INTRODUCTION 

In the last few years the study of physical phenomena 
characterized by a degree of self-organization M , has at- 
tracted a lot of interest. These models are usually cellular 
automata models defined on a discretized lattice, with a 
growth rule that can be either stochastic, when the in- 
homogeneities in the system change with a time scale 
smaller than the characteristic time scale of the dynam- 
ical evolution (noise), or deterministic with a quenched 
disorder which accounts for the effect of inhomogeneities 
inside a solid medium. Both kind of dynamical rules are 
characterized by an evolution towards an attractive fixed 
point in which scale free fluctuations in time and space 
are present |g| . 

The problem of kinetic roughening belongs to this class 
of models. It received recently an increasing interest in 
relation with non-equilibrium growth models |3j and in 
view of its practical applications: Chemical Vapor Depo- 
sition (CVD) H and electro-chemical deposition H are 
just two examples. 

In this perspective, one is interested in identifying 
the dynamic universality classes of kinetic roughening 
processes and several models has been defined starting 
from the models falling in the universality class of the 
Kardar-Parisi-Zhang (KPZ) equation H. This equation 
describes the properties of an interface h(x, t) driven by a 
stochastic noise and gives a roughness exponent x = 0.5. 
Other models are more suitable to describe the propaga- 
tion of interfaces in random media, i.e. with a quenched 
disorder. These models are driven by an extremal dy- 
namics. In this class fall the so-called Sneppen model J7|] 
(in |7[| referred as model B) and the pinning model by 
directed percolation || , which predict a roughness expo- 
nent equal to x = 0.63. These models produce self-affine 
surfaces. Recently, a model has been introduced to de- 
scribe some etching experiments, which leads to the for- 
mation of self-similar (fractal) structures, and which has 
been shown to fall in the percolation universality class 



||. Many experiments on surface roughening [|l0[— |l2|| , 

however, as well as experiments on chemical etching |13j ] 
produce self-affine surfaces instead of self-similar ones. In 
this paper, we will focus on kinetic roughening phenom- 
ena leading to self-affine (x < 1) surfaces. 

We recall that the roughness exponent is defined by the 
ensemble averaged width of the interface as W(l, t) =< 
(h{x,t)- < h(x,t) >) 2 >!/ 2 - l^f(t^ z /l) where z is the 
so-called dynamical exponent, the angular brackets de- 
note the average over all segments of the interface of 
length I and over all different realizations. / is a scaling 
function such that f(y) ~ y x for y << I and f{y) = cost. 
for y » 1. The exponent [3 = x/ z describes the tran- 
sient roughening, during which the surface evolves from 
the initial condition toward the final self-affine structure. 

In spite of the strong universality exhibited by the 
KPZ and Sneppen models, in many experimental stud- 
ies one measures values of x which are above the ones 
predicted by both the KPZ and Sneppen universality 
classes. To give some examples, we remember the ex- 
perimental studies of paper burning for which one gets 
X = 0.70 ± 0.03 pi, or the propagation of a forced fluid 
front in a porous medium, which exhibits a roughness 
exponent X = 0.73 ± 0.03 Q and x = 0.88 ± 0.08 |§. 

In this paper we propose a generalized model for ki- 
netic roughening characterized by anisotropic growth 
rules and, as a consequence, separated time scales for 
the dynamics. The existence of this time-scale separa- 
tion induces a cross-over effect in the roughness proper- 
ties, which could erroneously appear as a genuine non- 
universal critical behavior, and could give an explanation 
for the above cited scattered experimental results. 

Some results presented here have been already briefly 
reported on in a letter |jI3| . In this long paper we give 
a detailed, complete description of our previous work. 
Moreover, we present a set of new numerical results, 
which allow us to reach different and better founded con- 
clusions with respect to O] . 



The idea underlying the model is that some experi- 
mental parameters can introduce a characteristic scale 
in the system, separating different scaling behaviors. In 
particular we consider a model which includes explicitly 
an anisotropy factor, say a growth rule dependent on the 
local environment of the growing site. The model thus 
presents a complex interplay between a global equilib- 
rium and the conditions of a local dynamics. This choice 
is motivated by the observation of roughening phenom- 
ena occurring in etching processes which represent an im- 
portant tool either in academic research or in device tech- 
nology. Their importance is related to the preparation of 
single-crystal samples of desired dimensions, shapes and 
orientations. Etching is usually applied to obtain desired 
mesas and grooves in semiconductors wafers and multi- 
layers Q. 

In the same field, although in a different context, 
etching processes are used to produce textured optical 
sheets, which allow to exploit the light trapping by to- 
tal internal reflection to increase the effective absorp- 
tion in the indirect-gap semiconductors crystalline sili- 
con. Light trapping, originally suggested to increase the 
response speed of silicon photo-diodes while maintaining 
high quantum efficiency in the near-infrared, was later 
indicated as an important benefit for solar cells |L4J . 

The general suffix etching indicates the ensemble of 
operations which involve the removal of materials by ex- 
pending energy either by mechanical, thermal or chemical 
means. In ref. Jl3| the authors focused their attention on 
chemical etching processes as a reference point to formu- 
late the model. One of the most important properties of 
these processes is represented by the intrinsic anisotropy 
p5| , ^6| of the etch rates. For instance in samples of crys- 
talline silicon etched in solutions of aqueous potassium- 
hydroxide (K-OH) with isopropil alcohol (IPA), depend- 
ing on the concentration of the etchant and the temper- 
ature, the (111) direction etches slower than the others 
by a factor which can be of order 100 or more [O. The 
degree of anisotropy affects the properties of the surface, 
which turns out to be rough with an apparently non- 
universal roughness exponent. 

Although the definition of the model is very general 
we will briefly consider the chemistry of the etching pro- 
cess in order to exhibit a physical framework that allows 
to understand the meaning of the definitions and their 
interpretation. 

The disorder in the etching process is related to the 
impurities in the lattice. Such impurities, e.g. vacant 
atoms, reduce the binding energy of atoms nearby the 
vacancy. By assigning to each site (atom) of our lattice 
a random number Xi we assume that a distribution of 
vacancies, or other kinds of impurities, is present in the 
system, and this induces fluctuations in the binding en- 
ergy of atoms due to this disorder. If we assume to be in 
a condition of slow dynamics, that is to say the driving 
field (which in our case in represented by the concentra- 
tion of etchant) tends to zero jlq ], we can look at the 
etching as an extremal process, where the etchant dis- 



solves the atom with the smallest binding energy. This 
is correct for low etchant concentrations and corresponds 
actually to the situation experimentally more interesting, 
in which rough surfaces are produced. 

In order to reproduce the experimental conditions 
(type of etchant, concentration, temperature) a micro- 
scopic model for the physical process should contain some 
tunable parameters (at least one). In our model, the 
anisotropy is introduced by a phenomenological tunable 
parameter, p, which distinguishes sites with a different 
local environment. 

The introduction of the parameter p defines a char- 
acteristic scale in the problem. As a result the critical 
properties of the model are characterized by a continuous 
crossover between two universality classes corresponding 
to the roughness exponents x = 1 an d X — 0-63 (Sneppen 
models A and B respectively 0). In particular one can 
define a parameter r = j?— which measures the time- 
scale separation between the dynamics of the different 
classes of sites. For lengths I > I* <x - ~ p one observes 
a behavior characteristic of the Sneppen model B uni- 
versality class (x = 0.63), while for / < I* ex ' ~ p one 
observes a behavior characteristic of the Sneppen model 
A universality class. The existence of this crossover is dif- 
ficult to detect directly on the plot for the scaling of the 
surface width W 2 (l) (especially for finite sets of data) and 
it becomes evident looking at the power spectra. This 
explain why large-scale experiments could give the im- 
pression of non-universality in the critical properties of 
rough surfaces. 

Let us look at the meaning of p in the case of etching. If 
we represent the crystalline lattice of the silicon on a two- 
dimensional plane we can imagine a square lattice where 
(see Fig. ([!])) the atoms can be found in each of the four 
positions marked in figure by the letters a — d. The four 
positions correspond to different oxidation states: from 
the situation (c) (oxidation number 0) which occurs only 
in the bulk, to the situation (d) (oxidation number —3). 
Note that all the surface atoms are passivated by hydro- 
gen atoms. The atoms in the positions (a) and (6), corre- 
sponding respectively to the oxidation numbers —2 and 
— 1 (two and one heteropolar bonds, i.e. Si-H bonds), 
play an important role in explaining, at least from a 
heuristic point of view, the origin of the anisotropy in 
the etched rates |16J . The parameter p quantifies the ra- 
tio of the etch rates between the sites in the positions (a) 
and (b). The basic idea is that in the (111)— plane of sili- 
con, there is only one heteropolar bond per silicon atom. 
Therefore there are three bonds to break for dissolution, 
while other planes (except the (110)) have more than one 
heteropolar bonds and accordingly a smaller number of 
bonds must be broken. 

The paper is organized as follows. In section II we 
describe in detail the model and the setup of numerical 
simulations. In section III we present and discuss the 
numerical results in relation with the Sneppen model A 
and B universality classes. In section IV a discussion of 



the results and some conclusions are drawn together with 
a planning of future researches. 



II. THE MODEL 

We give now a detailed definition of the model. The 
model is defined on a square 2D lattice tilted at 45° 
(see Fig. [y). We consider a 1 + 1 dimensional inter- 
face h(x) = h(x,t) with x — 1,2,...,L, where L is the 
linear extension of the interface in the x direction. The 
initial condition for the dynamical evolution of this in- 
terface is given by: h(2x, 0) = 1 Vx G [l,L/2] and 
h(2x - 1,0) = Wx E [l,L/2] , in order to have both 
classes of variables (lattice planes) participating to the 
dynamics from the beginning, but different initial con- 
ditions do not change the properties of the model. The 
interface, which satisfies locally the conditions Il9| 



\h(x,t) + 1 - h(x - l,t)\ < 1, 
\h(x,t) + 1 - h(x + l,t)\ < 1, 



(1) 



contains two classes of random variables that correspond 
to two separate classes of sites. The sites (M) for which 
it holds V 2 /i > (called minimum sites) which are, mi- 
croscopically, the atoms with two heteropolar bonds, and 
the sites on a slope (slope (S) sites) for which one has 
V 2 /i = 0. These last sites correspond microscopically 
to atoms with one heteropolar bond. To each class of 
sites is assigned a class of Gaussian distributed uncor- 
rected random variables which mimic the disorder, and 
represents physically, for the case of etching, the binding 
energy of atoms: 

, , J [0 : 0.5] if x is such that V 2 h = . . 
r/{x, h) e | jj Q 5 . ^ if x ig guch that V 2 h > {<&) 

The sites with V 2 h < 0, for which all the chemical bonds 
are homopolar, i.e. Si-Si, do not take part to the dy- 
namics and they have assigned a zero value of the ran- 
dom variable. Periodic boundary conditions are assumed 
along the x direction. 

The system evolves by updating the site i* with the 
largest r.v. in one of the two classes of sites chosen, at 
its turn, with a probability p. One thus updates with 
probability p a site (S) and with probability 1 — p a site 
(M) according to the rules (see Fig. (0)): 

(1) h(i*,t+l) = h(i*,t)+2, r)(i*,h(i*,t+l),t+l) = 0; 

(2) Updating all the sites necessary to make satisfied 
the conditions [j] (this phase is assumed to be in- 
stantaneous with respect to extremal dynamics); 

(3) Updating of the random variables for the sites 
which changed their class of belonging. In partic- 
ular T)(x,h,t + 1) = 1/2 * RAN if r](x,h,t) = 
and T](x, h, t + 1) = T)(x, h, t) + 1/2 if t](x, h, t) ^ 0, 
where RAN is a random value between and 1; 



(4) Updating of the random variables of the sites which 
have changed their height but which did not change 
their class of belonging, (sites S only): r](x,h,t + 
1) = 1/2* RAN. 

The parameter p can vary in the range [0 : 1/2]. If we 
define t$ as the characteristic time scale for S variables 
and £m the characteristic time scale for M variables one 
has: 



ts 



P 



l-p 



(3) 



The growth of the interface, in the etching process, rep- 
resent the invasion of the etchants into the silicon wafer. 
^From this point of view the updating of the sites (5) 
mimics the etching of the (111) planes and the updating 
of the (M ) sites the etching in the (100) direction. For 
p = 1/2 all the sites which take part to the dynamics are 
equivalent and there is no anisotropy (r = 1), whereas 
the case p = corresponds to the maximal anisotropy in 
which f(iii)/f(ioo) = r = 0, where V(m) and W(ioo) are 
the etch rates in the corresponding directions. 

Our model can be viewed as a variation of the Sneppen 
model for quenched surface growth, where two important 
elements are added: 1) The anisotropy in the distribu- 
tion of the quenched random field, depending on the local 
characteristics of the growing surface: 2) a time scale sep- 
aration for the dynamical evolution of the two classes of 
variables (S sites and M sites), which is tuned by the 
parameter p. 



III. NUMERICAL RESULTS 



We have studied this cellular automata by numerical 
simulations in order to analyze its dynamical roughening 
properties. The sizes we have chosen for the numerical 
simulations range from L — 2048 to L — 8192. For each 
value of p(we have considered p = 0.0, 0.02, 0.2, 0.5) , 10 2 
simulations lasting 10 7 time steps have been performed 
and we have computed the growth exponent /3, which 
rules the time evolution of the width W(t) of the surface 
(W(t) ~ t@) before the stationary state is reached, and 
the roughness exponent \, which gives the scaling of the 
width of the surface (W(l) ~ l x ), in the stationary state. 
The stationary state is called self-organized in that it is 
reached spontaneously by the system independently of 
the initial conditions. This self-organization is confirmed 
by an analysis of the temporal evolution of the distri- 
bution of quenched variables (the histogram <& t (j))). To 
characterize temporal correlations in the dynamics and 
check that the asymptotic state is critical, we studied the 
distribution of the avalanches in the asymptotic state. As 
an independent check about the universality of the rough- 
ness properties of the model, we have studied numerically 
the power spectrum S(k) of the height profile. 



To ensure that the system is in the stationary state, we 
studied the behavior of the n-th moments (for n = 3,4, 5) 
of h(x,t) normalized to second momentum: 



m n (t) 



(Ei(h(i,t)-h(t)) n ) 

«Ei(fc(»,t)-ft(t)) a » J 



(4) 



where h(t) is the mean surface height at time t. We get 
that, after a transient, all the odd moments vanish and 
the even ones tend to constant values (see Fig.s |^-||). In 
particular the condition for the skewness 777,3 = (Fig. 
|3j), which characterizes the stationary critical state J/J, is 
realized after about 10 3 time steps per site, independent 
of the value of p. These results imply that the higher 
moments of the variable h(x, t) — (h) scale in a trivial 
way (they are powers of the second moment), and after 
the transient the probability distribution of the variable 
h(x,t), which can be viewed as a random variable, is 
Gaussian. The amplitude of the normalized even mo- 
ments m2m in the asymptotic stationary state, charac- 
terizes the roughness properties of the interface. 

In Table B we report the measured values for the dy- 
namical exponents (3 which turns out to be independent 
of p. Fig.s (p|||) show the scaling behavior of W 2 (l) 
for different values of p. For p = 0.0 (i.e. maximal 
anisotropy) the measured values of \ are affected by a 
finite-size effect and they tend, in the limit L — » oo, to 
the value x = 1.0 found for Sneppen model A j^J (Fig. ||). 
In this case the surface is composed by very big pyramids 
(Fig. |l0| ). On the other hand for p — 0.5 one recovers 
the universality class of Sneppen model B with \ = 0-63 
(Fig. |]). For all the other values of p between and 0.5 
(p = 0.02 in Fig. | and p = 0.2 in Fig. ||), trying to 
perform fits away from the saturation regions, one would 
be tempted to invoke the existence of a non-universal be- 
havior ruled by the parameter p. A careful observation 
puts in evidence that the curves for W 2 (l) seem to ex- 
hibit a crossover between the Sneppen models A and B 
universality classes. On the basis of p one can define a 
characteristic length /* oc ' ~ p above which one could 
see the \ — 0.63 behavior and below which the X = 1 
behavior. We shall come back on these considerations 
later on when we shall discuss the power spectra. 

We have also studied the time evolution of the distri- 
bution of random variables on the invading interface (the 
histogram <J>i (77) , where 77 is a generic value for rj(x,h)), 
which is of great importance for models with extremal dy- 
namics. The results of the simulation are shown in Fig.s 
(O-hJ). One can see that $4(77) self-organizes, for p =/= 0, 
after about 10 time steps, into a distribution that is the 
superposition of two theta functions, one for each class 
of variables, each one characterized by a critical thresh- 
old Pc(p) depending on the parameter p. The meaning 
of these thresholds is that only S variables larger that 
7?f and M variables larger than r^ 1 can grow M. For 
p = 0, instead, the histogram has no self-organized crit- 
ical state (Fig. |lTJ) . Looking carefully at Fig. ( pi] ) we 
can see that, while in the initial transient there arc a few 



S sites, in the asymptotic state most sites are S sites. In 
fact nearly all variables larger than 0.5 (the M variables) 
are disappeared. This observation agrees with the actual 
structure of the surface, which is composed by very big 
pyramids (Fig. n0|), with a roughness exponent x — 1- 
This picture is confirmed by the acceptance profile a(rf), 
which is shown in Fig. n5l As in the Bak and Sneppen 
model, the acceptation profile (that is to say the distri- 
bution of the values of all updated quenched variables up 
to the actual time) exhibits correlation properties (it is 
not flat), reflecting temporal correlations in the dynam- 
ics. But, while the acceptation profile for S variables is 
quite similar to that of the BS model, going to zero lin- 
early at p c , the acceptation profile for M variables has 
a more complicated behavior. This difference originates 
maybe from the fact that S variables (77 G [0.5,1]) can 
turn into M (77 e [0,0.5]) variables during the dynamics 
of the system, while M variables cannot become S vari- 
ables. Moreover, S variables can have developed correla- 
tions before the transition to M variable and this affects 
the shape of a(rj) for 77 s [0.5, 1]. This might account for 
the linear part of the acceptation profile of M variables, 
around 77 = 1, but the non linear part is more puzzling. 
The coupling between S variables and 777, variables could 
play a role in this behavior, too, but at the moment we 
have no clear explanation of it. ^From the a(j]) we can 
get a good estimation of the critical thresholds 77;? and 
77^ for different values of p (see Table [V). 

The stationary state is characterized by a constant ra- 
tio between S sites and M sites, that is to say the evo- 
lution equation for the densities ps and pu of sites S 
and M respectively, have an attractive fixed point in the 
stationary state (see Fig. fijla-d), with the asymptotic 
values of ps, Pm depending on the parameter p. One in- 
teresting observation is that, even in the case p = 0, that 
is to say only M sites can be selected by the extremal 
dynamical rule, there is a stationary state for the system 
with ps 7^ 0. This is due to the particular geometry of 
the lattice, for which the growth of an M site implies the 
creation or annihilation of some S sites. In other words, 
there cannot be surfaces without slopes (S sites). 

The roughness exponent accounts for scale free spatial 
fluctuations in the interface profile. In order to char- 
acterize the eventual scale free fluctuations in the dy- 
namical evolution of the system at its asymptotic crit- 
ical state, that is to say long range temporal correla- 
tions, we have studied the avalanche distribution. An 
avalanche is defined as a sequence of causally connected 
elementary growth events. For the class of models with 
quenched disorder and an extremal dynamics to which 
our model belongs, the initiator of a critical, scale in- 
variant, avalanche is identified in the critical state by 
a site with quenched variable r)^(p) or rjf{p) (respec- 
tively for an M initiator and for an S initiator). The 
values of 77^ and 77;? for different values of p can be ob- 
tained by the asymptotic histogram distributions shown 
in Fig.s (111- 14). In our case there are two classes of vari- 
ables, the S and M sites, and two possible initiators for 



an avalanche. We call the avalanches that start with an 
S site, S-avalanches, and the avalanches that start with 
an M site, M-avalanches. An avalanche lasts when a 
variable which has been updated before the growth of 
the initiator is selected by the extremal dynamics. The 
statistics of off-critical avalanches has been shown to have 
the form |@: 



P*(s;r,) = s- rx f x (\r,-r,t\s (rx ) 



(5) 



where X = S, M , and rj is the initiator of an X-avalanche. 
This distribution becomes a pure power law for r\ = r^ . 
In the limit t — > oo the system self-organizes into the 
critical state r\ = r^, and the (normalized) avalanche 
size distribution becomes: 



P 



x 



(rf) = 



e: 



(6) 



We have performed a set of about 10 3 realizations of size 
L = 8192, lasting each one 2 x 10 6 time steps, and col- 
lected the statistics of S and M avalanches over the last 
10 6 time steps, for p = 0.02,0.2,0.5. These simulations 
required about 2 months of CPU time on our comput- 
ers (a network of DEC alpha machines with clocks going 
from 266MHz to 500MHz), and are at the best of our 
computation possibilities. To reduce numerical problems 
connected with the approximation on r]^f , rjf , we used 
an alternative definition of critical avalanches in models 
with extremal dynamics, which resides on the causal re- 
lation between updated sites inside an avalanche (for de- 
tails on the definition of critical avalanches see [glj f23[ ) . 
The results are shown in Figs. (|17|), (jig). Even after this 
big computational effort, our numerical results are still 
a bit noisy. In particular the statistic of S avalanches 
for p = 0.02 is really poor. This is due to the fact that 
for small p values most of sites selected by the dynam- 
ics are M sites. Consequently, it is difficult to observe a 
quite clear power law behavior for both the S-avalanchcs 
and M-avalanches distributions. We point out that the 
presence of long range temporal correlations is not neces- 
sary for the model to have self-similar or self-affine spatial 
properties, as already observed in a different context P6| . 
In order to better establish the critical properties of 
our model we have measured the power spectrum S(k) 
of the equilibrium surface. The model studied here is a 
discrctized cellular automaton which can be thought as 
a modified version of the Sneppen model for quenched 
interface growth. The Sneppen model has been shown to 
be, at least in 1 + 1 dimensions, in the same universal- 
ity class of the continuous Kardar Parisi Zhang equation 
with quenched noise (QKPZ) p4| . It is natural, but not 
necessarily true, to suppose that for our model, too, it 
is possible to find a formulation as a continuous growth 
equation. Given a general growth equation for h(x,t) like: 



dh(x,t) 

at 



= A[h(x,t)] +i(x,t) 



(7) 



where A[...] is an operator acting on h(x,t) and j(x, t) 
is an uncorrelated quenched noise (the "temporal" direc- 
tion corresponds to the growth direction of the surface) , 
if the operator A[...] is linear and local, the equation can 
be Fourier transformed into 

iujh(k,w) = A(k)h(k,u>) + j(k,u>), (8) 

and by introducing the propagator G(k,u), 



h(k,tu) — G(k,oj)j(k,u>). 



(9) 



where G(k,u) = [iw — A(k)]^ 1 . 

The propagator G(k, u>) of eq. is related to the power 
spectrum S(k) of the interface in the asymptotic state. 
The power spectrum is so defined: 



S(k) = {FT[h(x)h(x')]) = (\h(k)\ 2 ) 



(10) 



where FT[...] is the Fourier transform operator, the av- 
erage is over different realizations of the noise, h(x) — 
h(x, t = oo), h(k) is the Fourier transform of h(x). Eq. 
[lOl is valid is the case the noise is uncorrelated in space 
and time, which is the case of our model. The relation 
between G(k,u>) and S(k) is the following ( 0): 



G(k,w = 0) 2 = S(k) 



(11) 



Equation (111]) tells us that the power spectrum of the in- 
terface can give informations on fc-dependent part of the 
propagator G(k,cu = 0) and consequently on the struc- 
ture of the operator A(k) in Eq. ||. For self-affine sur- 
faces, the power spectrum follows a power law scaling 



S(k) 



-2<5 



(12) 



where 8 is related to the global roughness exponent Xgiob 
through the scaling relation 28 = 2xgiob + 1 fcj • 

Fig.s dl9-E0h report the behavior of the power spec- 
trum S(k) of the interface profile in the critical state for 
system sizes L = 2048, 8192 and different values of p. 

For large values of k finite size effects connected to 
the discretized nature of the model become relevant 
and there is a deviation from the power law behavior. 
Away from this saturation effect it is evident in this 
case how S(k) is characterized by a clear cross-over be- 
tween two power law behaviors. In the low k region 
(k < k* (x p/(l/2 - p) = 1//*) S(k) scales with an 
exponent 8i ow close to 1 (actually 1.07(3)) for all val- 
ues of p. The QKPZ (Sneppen B model) universality 
class is characterized by a global roughness exponent 
(coinciding with the local exponent we have measured 
through W 2 {1)) Xgiob = X = 0.63, giving S = 1.13, 
quite near to what we find. For intermediate values of 
k (k > k* (x p/(l/2 — p) = l/l*) one gets an exponent 
8 m id quite near to the value 1.7 (we find values between 
1.7 and 1.86) which corresponds to the QEW universality 
class {xgiob ~ 1-2 [gj), independently of the value of p. 
We point out that the QEW model is super rough with a 



global roughness exponent Xgiob ~ 1.2, and a local expo- 
nent x = 1-0 (the one we measured through W 2 (l)). For 
p near 0.5 the intermediate k region is very small, and 
it is difficult to distinguish it from the saturation region. 
The same happens when one tries to fit the low k region 
for p close to zero. lip = 0.0 or = 0.5, no cross-over effect 
is observed, since these two values of p correspond to the 
"pure" QEW and KPZ universality classes, respectively. 
From these results we have a confirmation that the 
model does not exhibit non-universal critical properties. 
The apparent non-universal roughness exponent is the 
consequence of a cross-over effect, tuned by the param- 
eter p. The fact that this cross-over effect is difficult to 
observe when studying the scaling of the mean square 
width W 2 (l) of the sample, could explain the discrep- 
ancy between the experimental findings available in the 
literature. 



and the anisotropy effect which takes into account local 
constraints in the growth. 

It is worthwhile to stress how our model suggests the 
possibility of several analytical approach, from the treat- 
ment of the problem in terms of a continuous stochas- 
tic dynamical equation, to the single site mean-field ap- 
proach p7| , or to the application of a method recently 
proposed for dynamical models driven by an extremal 
dynamics |2l] |2^j2^|. Particularly promising, in this re- 
spect, is a recently proposed non-perturbative Renormal- 
ization Group approach [M which allows one to study 
self-affine problems. 

The author acknowledges financial support under the 
European network project FMRXCT980183. 



IV. CONCLUSIONS 

In this paper we have introduced a model for surface 
roughening whose main peculiarity is that of taking ex- 
plicitly into account the anisotropy of the growth pro- 
cess by means of a tunable phenomenological parame- 
ter p which introduces local, i.e. dependent on the local 
environment, dynamical rules in the growth. The sim- 
ple introduction of just one anisotropy parameter p is 
far from being able to capture all the characteristics of 
etching processes, and in general of surface roughening 
experiments. In etching experiments, for example, trans- 
port phenomena in the solution are likely to be impor- 
tant and both concentration and agitation have strong 
effects on transport. Nevertheless, our model captures 
at least some basic elements of the relationship between 
anisotropy and the apparent non-universality observed 
experimentally in etching processes. Moreover, the gen- 
eral requirement of a microscopic dynamical rule depend- 
ing on the local environment could be a key element in the 
apparently observed non-universality in kinetic roughen- 
ing phenomena. 

As a main outcome, the model exhibits a cross-over 
behavior in its critical properties. For each value of the 
anisotropy factor p the system reaches a critical station- 
ary state, with a characteristic length separating a KPZ- 
like (Sneppen B model) behavior from a QEW-like (Snep- 
pen A model) behavior. The cross-over from one scaling 
behavior to the other is tuned by the anisotropy parame- 
ter p. If one looks at the scaling of W 2 (l), the cross-over 
effect cannot be easily discovered, and the system seems 
to have a non- universal roughness exponent. A careful 
analysis of the power spectrum S(k), however, shows a 
clear cross-over effect. These results can probably help 
to explain the relevant discrepancies among experimen- 
tal results |1Q-|F4]. We believe that this behavior is the 
outcome of the complex interplay between the global dy- 
namics which selects at each time step the weakest site 
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SILICON 
FIG. 1. Schematic representation of the crystalline silicon 

lattice as a square lattice. 
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FIG. 2. Schematic representation of interface dynamics. 
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FIG. 4. Time evolution (adimensional units) of the mo- 
ment m,4 (adimensional units) of the growing interface, nor- 
malized by the second moment, for p = 0.0, 0.02, 0.5. One sees 
that asymptotically m.4 tends to different constant values for 
the different values of p. 
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FIG. 5. Time evolution (adimensional units) of the mo- 
ment m,5 (adimensional units) of the growing interface, nor- 
malized by the second moment, for p = 0.0, 0.02, 0.5. One 
sees that asymptotically 7715 vanishes for all values of p. 
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FIG. 3. Time evolution (adimensional units) of the mo- 
ment m3 (adimensional units) of the growing interface, nor- 
malized by the second moment, for p = 0.0,0.02,0.5 (skew- 
ness). One sees that asymptotically 7713 vanishes for all values 
of p. 
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FIG. 6. W 2 (l) vs. I (both W 2 (l) and I are expressed in 
adimensional units) for p — 0.0 and different system sizes. 
The lower fit (dot-dashed line) corresponds to a size L = 512, 
giving and exponent x = 0.88(2), while the upper fit (dashed 
line) corresponds to a size L = 8192 and gives an exponent 
X = 0.96(2). In this case we expect that the exponent con- 
verges to x = 1 m the limit L — > co. 
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FIG. 7. W 2 (l) vs. I (both W 2 (l) and I are expressed in 
adimensional units) for p — 0.02 and different system sizes. 
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FIG. 8. W 2 (l) vs. I (both W 2 (l) and I are expressed in 
adimensional units) for p — 0.2 and different system sizes. 




FIG. 9. W 2 (l) vs. I (both W 2 (l) and I are expressed in 
adimensional units) for p = 0.5 and different system sizes. 




FIG. 10. A realization of the growing surface for p = 0.0 
(horizontal adimensional position on the x axis). The sur- 
face is composed by very big pyramids, thus with a strong 
prevalence of S sites. 
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FIG. 11. Histogram $(??) (r/ is an adimensional number) 
of quenched variables, at different times t (adimensional com- 
puter units), for p — 0.0. Asymptotically, all (most of the) M 
variables are eliminated. 
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FIG. 12. Histogram $(ri) (r) is an adimensional number) of 
quenched variables, at different times t, for p — 0.02. 
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FIG. 13. Histogram $(77) of quenched variables, at different 
times t, for p = 0.2. 
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FIG. 14. Histogram $(??) (r] is an adimensional number) of 
quenched variables, at different times t, for p — 0.5. 
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FIG. 16. Time evolution (£ is expressed in adimensional 
computer time units) of the densities ps and pm of sites S 
and M respectively, for p = 0.0 (a), p = 0.02 (b), p = 0.2 (c), 
andp = 0.5 (d). 




FIG. 17. Binned S-avalanches distribution for different val- 
ues of p. 




FIG. 15. Asymptotic acceptation (not normalized) profile 
a (v) (v is an adimensional number) for p = 0.02 (circles), 
p — 0.2 (squares) and p — 0.5 (triangles). 
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FIG. 18. Binned M- avalanches distribution for different 
values of p. 




FIG. 19. Power spectrum S(k) (all the quantities are ex- 
pressed in adimensional units of the computer simulation) of 
our model for p — 0.0, 0.02, 0.2, 0.5 (values referring to, re- 
spectively, the plots from bottom to top) and L = 2048. As a 
guide for the eye, we report the scaling law for KPZ (dotted 
line) and QEW (dashed line) universality classes. 




FIG. 20. Power spectrum S(k) (all quantitities are ex- 
pressed in adimensional units of the computer simulation) of 
our model for p — 0.0, 0.02, 0.2, 0.5 (values referring to, re- 
spectively, the plots from bottom to top) and L = 8192. As a 
guide for the eye, we report the scaling law for KPZ (dotted 
line) and QEW (dashed line) universality classes. 
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TABLE I. Values of the dynamical exponent j3 in our 
model for different values of the anisotropy parameter p. 
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0.02 
0.2 
0.5 


0.47(1) 
0.41(1) 
0.35(1) 
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0.54(1) 
0.63(1) 



TABLE II. Critical thresholds r), 
M, for different values of p. 
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